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I. INTRODUCTION 

Phospholipid membranes and self-assembled bilayers 
of block copolymers have many common features orig- 
inating from the amphiphilic nature of the molecules. 
Phospholipids comprise of hydrophilic head group and 
two hydrophobic acyl chains and the driving force of self- 
assembly in membranes, as in case of block copolymers, 
is the hydrophobic effecti Thus, theoretical methods de- 
veloped for the self-assembly of block copolymers can be 
applied for the description of phospholipid membranes 2 . 

Theoretical description and computer simulation of 
phospholipid bilayers and biological membranes is a chal- 
lenging task for many decades 2 ^. Although phospho- 
lipids are relatively short chains in comparison to poly- 
mers and block copolymers, the configurational space of 
conformations and the number of interactions are ex- 
tremely large. In addition, the simulation of phospholipid 
bilayers and membranes is a problem involving many in- 
teracting molecules, thus full description of collective be- 
havior would require simultaneous simulation of a large 
number of molecules. That is why, full atomistic molec- 
ular dynamics models^—, describing membranes with 
chemical accuracy, are limited to short times and short 
lengths. Usually, atomistic models on times exceeding 
microseconds and lengths larger than nanometers are not 
practical^, especially when the equilibrium structures of 
large systems are considered. 

One of the strategy to overcome this problem is to 
group atoms into effective particles that interact via effec- 
tive potentials within so-called coarse-grained models^. 
Coarse-graining can reduce considerably the degrees of 
freedom and the configurational space such that larger 
systems and longer times can be reached. Monte- 
Carlo (MC) simulations of phospholipid membranes with 
coarse-grained models of phospholipids^^— are proved to 
capture the essential properties of self-assembly of phos- 
pholipids into membranes, their molecular structure and 
collective phenomena. 

However, if the objective of the research is the equi- 



librium properties and the molecular details of the self- 
assembled structures, the equilibration time may largely 
exceed microseconds and thus, may not be reached with 
MD simulations in a reasonable time. Alternative strat- 
egy is to reach directly the minimum of the free energy 
using methods based on mean field theories, thus, avoid- 
ing computationally expensive equilibration process used 
in MC and MD simulations. 

The method of the Single Chain Mean Field (SCMF) 
theory 2 ^— is particularly suitable for such purposes: it 
combines analytical theory of the mean field type with 
the conformational variability accessible in MC simula- 
tions. Thus, the accuracy of the method in describing 
the molecular details of equilibrium structures can com- 
pete with MC simulations while the speed of obtaining 
the results depends only on the speed of the numerical 
solution of equations, thus it can be much faster than 
MC simulations. 

The advantage of this method is the speed in local- 
ization of the free energy minimum and the equilibrium 
properties as well as the precision in the measurement 
of the equilibrium free energy which is not straightfor- 
wardly accessible with MD and MC simulations. The 
disadvantage is the mean-field nature of the solution that 
does not include fluctuations and inter-particle correla- 
tions. Nevertheless, the shortcomings of the approach 
can be compensated by the combination of the SCMF 
theory with MD simulations which can complement the 
method. Fast mean field method provides for equilib- 
rium structure with molecular details which, in turn, can 
be used as input initial configuration for MD simulations 
and provide the lacking dynamic information and fluctu- 
ations. 

The SCMF theory, originally developed for the mi- 
cellization problem of low-molecular surfactant o 20 i 21 , de- 
scribes a single molecule in the molecular fields. The 
multichain problem is reduced to a single chain in the 
external self-consistent field problem. The position, ori- 
entation and configuration of the molecule depend on the 
surrounding fields while the fields depend on the con- 
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figurations of the molecules, and finally the equilibrium 
fields are found self-consistently by numerical solution of 
a system of nonlinear equations. This method gives a de- 
tailed microscopic information on the configurations and 
averaged positions of the molecules, the optimal shape 
and structure of self-assembled structures, the distribu- 
tion of molecules in the aggregates, the critical micellar 
concentrations as well as the optimal aggregation num- 
ber and the size distributions of the micelles and the 
aggregate a 23 i 24 . This method is quite universal: it can 
be applied to solutions of linear or branched polymers, 
solutions of low-molecular weight surfactants and vari- 
ous additives, mixtures of various components and struc- 
tural and shape transitions. However, the limiting factor 
restraining the use of the SCMF theory is the computa- 
tional realization of a stable code that can efficiently solve 
self-consistent equations. There arc no standard pack- 
ages with universal implementation of the SCMF theory, 
which could make the use of the method accessible for 
a large number of people like Mesodyn packag o 26 i 27 for 
mesoscopic modeling of phase separation dynamics or an 
integrated simulation system for soft materials (OCTA) 
for multiscale modeling 2 ^—. Most of the previous works 
on SCMF theory^— use iterative methods or external 
mathematical libraries, which arc not optimized for a 
given problem and thus, computationally unstable and 
highly demanding in computer resources, especially in 
RAM memory. Thus, problems involving long molecules, 
large systems and geometries more than ID confront with 
the limit of available computer resources and thus, mak- 
ing the solution of the problem without special technical 
skills to be a quite difficult task. A modification of the 
method, Single Chain in Mean Field calculations 3 -! - — , 
avoids the necessity to solve self-consistent equations. 
The direct solution of the self-consistent equations is re- 
placed by MC equilibration of the chains in the quasi- 
instantaneous fields maintained at the most recent val- 
ues and updated after a predetermined number of MC 
moves. In practice, the MC equilibration would slow 
down the calculation with respect to direct solution of 
the equations, but it is still much faster than the direct 
MC method. 

It is noteworthy that SCMF method is similar in spirit 
to one of the first mean field model of the phospholipid 
membranes by Marcelj a 34 ' 35 constructed for modeling of 
the fluid - gel phase transition and based on phenomeno- 
logical potential of the Maier-Saupe type between the 
segments of the tails. The boundary between the solvent 
and the membrane is modeled as a planar surface with 
the phospholipid tails attached with a fixed grafting den- 
sity. Each tail interacts with the neighbors via mean field 
which is found sclf-consistcntly. 

In this work we report a computational tool provid- 
ing relatively fast and stable solution of the equations 
of the SCMF theory in different geometries and dif- 
ferent molecule structures. The SCMF theory is ap- 
plied to model phospholipid membranes. We show that 
three different coarse-grained models of phospholipids 



can adequately describe the equilibrium properties and 
the molecular structure of phospholipid membranes. 

The paper is organized as follows. Theoretical prin- 
ciples of the SCMF theory and main equations of the 
method are introduced in the section [TT] in the most gen- 
eral way. This theory is then applied to the simulation 
of phospholipid bilayers and three coarse-grained mod- 
els of phospholipid molecule are introduced in section 
Mil These models and the resulting equilibrium phos- 
pholipid bilayers are compared with experimental data 
and between each other in section IIVI Last section [V] 
summarizes the obtained results while the computational 
details which are necessary for the implementation of the 
method are described in the Appendix [AJ 

II. THEORY 

The SCMF theory is an example of the Self-Consistent 
Field (SCF) method where a single chain is described at 
the molecular level while the interactions between differ- 
ent chains are described through a mean molecular field 
which is found self-consistently. 

The conformations of a single chain are generated with 
the Rosenbluth algorithm 3 ^ or MC simulation a 13 ' 37 and 
the intra-molecular interactions are calculated exactly us- 
ing the model potential for interactions between the seg- 
ments of the chain (see Appendix IA II for details) . The 
probabilities of individual chain conformations depend on 
the mean molecular fields, while the values of the fields 
are calculated as the average properties of individual con- 
formations. The resulting equations for the mean fields 
and the probabilities of the conformations arc solved self- 
consistently with the original method described in the 
Appendix IA 31 The solution of these equations gives the 
equilibrium structures and the concentrations profiles of 
all components in the system as well as the most prob- 
able conformations of individual molecules. The power 
of this method is the speed in obtaining solutions (faster 
than MC and much faster than MD simulations) and the 
precise calculation of the energies. 

This method is quite universal and can be applied 
for the mixture of an arbitrary number of molecules of 
different types interacting with each other through the 
mean fields. The free energy of a system containing 
iVi, AT 2 , Nm molecules of types 1, 2, M can be writ- 
ten as a sum of three terms, the entropy, intra- and inter- 
molecular terms, all written in terms of kT, 

F = — (S) + (H intra ) + (H inter ) (1) 

where the angular brackets denote the average over the 
probability distribution function (pdf) p of the system, 
(...) = J . . . palp. This function contains all informa- 
tion about the equilibrium state of the whole system. If 
there are no strong correlations between interacting par- 
ticles, such as ion pair formation or bonds formation, we 
can neglect the correlations between the molecules and 
use the mean field approximation: the pdf p of many 
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component system factorizes into the product of single 
molecules pdfs p a (T a i)i 



M N a 



P~\ ~[p a {T al ) 



(2) 



a=l i=l 



where r Q j is the conformation of i-th molecule of type 
a. Such factorization into contributions of individual 
molecules allows us to derive a close set of equations for 

(r«ri). 

The entropy term in the expression ([T|) is written as 



(5) = -(lnpA) 



(3) 



where the factor A is the de Broglie length which has a 
quantum mechanics origin, 
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Since the constants A Q do not appear in the final ex- 
pressions, they can be treated as unknown normalization 
constants. Thus, after factorizing the pdf p (J5J), the en- 
tropy term (j^J) reads 
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(5) 



where the brackets on the right hand side denote the av- 
erage over the single molecule pdf p a . Similar arguments 
allow us to write the intra-molecular energy of the sys- 
tem as a sum of contributions of the single molecules of 
different types 
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Before writing a similar expression for the inter- 
molecular part, we assume that molecules of each type 
a comprise of subunits with different chemical structure, 
and thus different energy parameters. Subunits can rep- 
resent Kuhn segments in the chain or, in case of coarse- 
graining description, the groups of atoms or beads in a 
coarse-grained model. Thus, the interaction energy be- 
tween a molecule of type a in the conformation state T a 
and a molecule of type /? in the conformation state Tp, 
can be written as a sum over the types of beads a as 



(T a ,Tp)=Y^ J K(T ai r)c a p(Tp,r)dr 



(7) 



In this expression Cp(Tp,r) is the concentration of units 
of type a at the point r of a molecule of type (3 in the 
conformation state Tp. These units interact with the 
field u^(T a ,r) created by the molecules of type a. Thus, 



factorizing the pdf p ([2]) one can write the inter-molecular 
interaction free energy in the form 

M 

£ N a (N p -S a p)Y, / «{r))(c a p{r))dr 

a,p=l a •* 

(8) 

where 5 a p is the delta symbol. Thus, the interaction free 
energy is represented by the interaction of the average 
concentration of beads with the average fields at each 
point. 

The free energy is usually coupled with the incompress- 
ibility condition implying that the sum of the concentra- 
tions of all components in a solution is fixed. However, 
this condition implies the hard core repulsion between 
all beads in the system, and thus creates computational 
problems in converging the SCMF equations. To over- 
come this problem, we introduce the explicit incompress- 
ibility condition in every point r in the system, 
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Q = l 



where {(f> a (r)) is the average volume fraction occupied by 
a molecule of type a in the point r, while 4>q is the total 
volume fraction occupied by the molecules of all types. 

Combining all terms ([5]), ©, flSJ) together with the in- 
compressibility condition ([3]), the free energy of the sys- 
tem can be written as 
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where A(r) is a Lagrange multiplier and all the averages 
are taken over the single- molecule pdfs p a . Minimization 
of this functional with respect to p a gives the pdf of a 
single molecule of type a 
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exp \-H™ tra (r a )-J2(N0-S a p)x 



J u a a {Y a ,r)(c a p {r))dr + J \{r)4> a {T a ,r)dr\ (11) 

where Z a is a normalization constant which is found from 
the normalization condition J /5 Q (r Q )(ir a = 1. This con- 
stant has a meaning of a partition function of the system 
at equilibrium and — ln Z a is the total free energy of the 
system at equilibrium. 
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Once the average concentrations {c^(r)^ and the vol- 
ume fractions {4>p{r)) are known, this expression allows 
to calculate the probabilities of each conformation p a (T Q ) 
for all molecules in the system. In turn, if the probabil- 
ities /o a (r a ) are known, the average concentrations and 
volume fractions, being the the molecular fields in this 
problem are found from the self-consistency conditions, 



(<£(r)> = J c«(r a ,r) Pa (T a )dT Q 
((f> a (r)) = J <j) a (T a ,r)p a {r a )dT a 



(12) 



representing the averages over the probabilities of con- 
formations. 

The probability of each conformation can be written 
as p a (T a ) = -j- exp[-if e //(r Q )], where H e f f (T a ) is the 
effective Hamiltonian given by (fTTj) . which describes the 
system at equilibrium. The last term in the effective 
Hamiltonian corresponds to the steric repulsion of beads 
of all types. If one of the components is a one bead sol- 
vent, the only degree of freedom of the solvent molecules 
is the position in space, r. Hence, the volume fraction 
of the solvent, can be found from the incompressibility 
condition (El) 



4> s {r) = 0o - 
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In addition, the Lagrange multiplier A(r) can be ex- 
pressed through the concentration of the solvent. The 
pdf of the solvent p s (r), the concentration c s (r) and the 
volume fraction occupied by the solvent <p s (r) are related 
via the following expressions 
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where N s is the number and v s is the volume of the sol- 
vent molecule, while the molecular field 



cj) s (r,r') w v s S(r-r') 



(15) 



where 8{r—r') is the Dirac delta-function. Substitution of 
(fT4|) and (fT5"| into (jTTJ) gives the approximate expression 
for the Lagrange multiplier A(r) 

M . 

v s \(r)*ln<t> s (r)+Y,( N e- S ^H / <(r,r') {c%{r')) dr' 

P=l a J 

(16) 

where we have assumed that H l s ntra (r) = and omitted 
few constants, which will cancel out by the normalization 
of p a . 

It is noteworthy, that in our equations the fields 
((f> a (r)) and (c°(r)) formally are not related although 
they both correspond to the concentration of monomers. 
They are indeed related through the volume of the 



monomers only if the molecules are composed of non- 
overlapping beads, when the distance between the cen- 
ters is larger than the the diameter. However, if we want 
to conserve the possibility of describing the overlapping 
beads, the coefficient of proportionality between (<fi a (r)) 
and (c°(r)) is not known a priori and we keep them as 
independent variables. 

The equations ([TT ]) .([T2 ]t .(fT5 ]) and ||T5J) form a closed set 
of non- linear equations of the SCMF theory. The solution 
of these equations gives the equilibrium structures, the 
self-consistent molecular fields such as the concentration 
profiles of the beads of each type and the distribution of 
the solvent, and the probabilities of each conformation 
of the molecules in the fields and the accurate measure 
of equilibrium free energies. The implementation of the 
computational method as well as the technical details 
related to the solution of the equations are present in 
Appendix [A] 



III. APPLICATION TO PHOSPHOLIPID 

MEMBRANES 

The method of the SCMF theory described so far is 
applied to model the equilibrium structures and prop- 
erties of phospholipid membranes. The SCMF theory 
can provide the detailed information on the microscopic 
structure of the phospholipid layer such as concentra- 
tion profiles of all groups of atoms of the phospholipid 
molecule, the thickness of the membrane, the average 
area per phospholipid head group and the mechanical 
properties such as compressibility of the membrane and 
the surface tension. All these parameters can be mea- 
sured in the experiment^— and thus, allowing the direct 
comparison with experimental data. 

We have developed the C++ code which runs in par- 
allel using OpenMP shared memory platform on the 32- 
core AMD machines (see Appendix I A 3[) . To simulate 
a flat phospholipid layer we use one dimensional geom- 
etry and discretize the space into parallel cells with the 
same value of the average fields. The periodic boundary 
condition is applied together with the assumption of the 
symmetry of the layer with respect to the central plane. 

Despite the complexity of the structure of the phospho- 
lipid membranes, the driving force for the self-assembly 
of phospholipids into bilayers is the amphiphile nature of 
phospholipids. That is why, wc only consider the units 
of two types, hydrophilic, H, and hydrophobic, T. The 
units of both types interacts with each other through the 
square well potentials. We assume that all units, T and 
H, and the solvent molecules S, have the same size and 
the same interaction range, while the interaction ener- 
gies arc different. We consider three types of interac- 
tions, interactions between hydrophobic units, T-T, in- 
teractions between hydrophobic units and solvent, T-S 
and hydrophilic units with solvent, H-S. In the following 
we will show that such assumptions are justified for quan- 
titative study of phospholipid layers and, in principle, 
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FIG. 1. (a) Chemical structure of the DMPC phospholipid 
molecule. 44-beads (b), 10-beads (c) and 3-beads (d) mod- 
els of phospholipid molecules used in the present calculations. 
Blue beads correspond to hydrophilic monomers, red beads 
correspond to hydrocarbon monomers in the tails and orange 
beads correspond to hydrophobic monomers in the heads. 
The interaction parameters of red and orange beads are the 
same. 



TABLE I. Parameters used for the simulation of phospholipid 
bilayer for three models of DMPC phospholipid molecule. 





44B 


10B 


3B 


Units radius (A) 


1.90 


2.50 


4.05 


Interaction range (A) 


5.70 


7.50 


12.15 


T-T contact energy (kT) 


-0.40 


-1.50 


-2.10 


T-S contact energy (kT) 


0.00 


0.00 


0.00 


H-S contact energy (kT) 


-0.10 


-0.20 


-0.15 


Bond length (A) 


1.5 


5.0 


10.0 


Occupied volume fraction <f>o 


0.700 


0.675 


0.675 


Sampling (number of configurations) 


3 x 10 6 


10 6 


10 6 


Simulation box size (A) 


120.0 x 120.0 x 62.7 



cohere with approximations used in successful models of 
phospholipid membranes (see e.g. Ref. Il~9l ). However, 
if necessary, more types of beads and more complicated 
potentials can be implemented. 

The chemical structure of the DMPC molecule is de- 
picted in Figure QJi. We present three models of DMPC 
phospholipid molecule at a different level of coarse- 
graining. The first, most detailed model, represents the 
DMPC phospholipid as a two tails molecule of 44 beads 
(Figure [Tp). A carbon group CH 2 of the molecule is 
assigned to one T bead. A phosphate group, P0 4 , is 
represented by 5 H overlapping beads, placed at the ver- 
texes and center of a tetrahedron. The choline group 
NC 4 H 11 is represented in a similar way: one H bead in 
the center of a tetrahedron is surrounded by 4 T beads 
placed in the corners. The COO groups are represented 
by 2 H beads. The angles between most of the bonds are 
fixed at value 120°. The torsion angles the sequence of 
groups CH 2 — CH 2 — CH 2 are allowed to have only three 



fixed values, 0°, 120° and 240°, which corresponds to 
cis- and trans- conformations of the groups. The param- 
eters of the model (Tabic HJ were adjusted with a series of 
fast simulations (conformational sampling size ~ 30 — 100 
thousands and 40 layers in the simulation box) while the 
accurate data were obtained with several millions of con- 
formations and hundred layers in the box. 

Two others, more coarse-grained and less detailed, 
models are depicted in Figures [Tp andQJl. One of them 
represents two tails phospholipid with 10 beads and an- 
other is simply 3-beads freely joined together. In contrast 
to 44-beads model, the conformational space of these 
models is much more restricted and we need less sam- 
pling to produce accurate results. Thus, the simplest 
3-beads model is the fastest and demanding less com- 
puter resources phospholipid model which can simulate 
the self-assembly of phospholipids into the bilayers with 
realistic properties of DMPC phospholipid membranes. 

It is important to note, that we simulate a finite size 
box with the fixed volume V and the fixed number of 
phospholipid molecules Ni with the free energy F given 
by the expression ()10l) . However, wc can extrapolate our 
data to much larger system with larger number of phos- 
pholipids Nf and larger number of solvent molecules. 
The description of a larger system is possible, if we do 
not take into account any energy contributions related to 
the perpendicular undulations and bending of the mem- 
brane. In this case the simulation box represent a self- 
similar part of a larger system. The free energy of the 
larger system F* is in a simple relation with the free 
energy per lipid / 

F* const + N*f (17) 

which, in turn, is related to the free energy of the simu- 
lation box minus the entropy of the solvent 

fNl=F - V *lfrfo. (is) 

v s v s 

where 4>q/v s is the concentration of the pure solvent. 

Thus, the calculation of free energy of the box F with 
different number of molecules in the simulation box gives 
the free energy of the membrane per lipid /. The min- 
imum of the free energy per lipid / corresponds to the 
equilibrium density of the membrane, thus the second 
derivative of this energy with respect to the area per 
lipid in the minimum gives the compressibility modulus 
of the layer 

K = 2A eq f'(A eq ) (19) 
representing a measure of the rigidity of the membrane. 

IV. RESULTS AND DISCUSSION 

The numerical implementation of the method in the 
form of two structurally independent modules: the gen- 



FIG. 2. Typical mean field "snapshots" (set of most-probable conformations of the lipids) of the phospholipid membrane within 
the 44-beads model. 



eration of a single molecule conformations and the solu- 
tion of equations (see Appendix [5} , permits us to con- 
sider several models of the molecule at the same time. 
Thus, we can compare our three models of phospholipids 
and test their performance in reproducing the thermody- 
namic properties of phospholipid membranes. 

The most detailed 44-beads model (Figure [TJ3) results 
in equilibrium structures of membranes where the phos- 
pholipids are assembled into fiat bilayers in fluid phase. 
Since the output of the method is the probabilities of each 
molecular conformation in the self-consistent fields, the 
most probable conformations may be used for the visu- 
alization of the resulting structures which correspond to 
the solutions of the equations. Typical mean field "snap- 
shots" of the equilibrium bilayer structures obtained with 
the 44-beads model are present at Figure [5J Although 
the picture visually resemble instant snapshots of MC 
simulations, this is not a representation of interacting 
molecules at a given moment of time, but these are the 
most probable conformations of a single chain in the fields 
corresponding to the equilibrium solution. Such "snap- 
shots" are helpful for visualization purposes in order to 
distinguish between different molecular structures. 

The equations of the SCMF theory may result in sev- 
eral sets of solutions corresponding to various equilibrium 
and even metastablc structures. To this end, the numeri- 
cal method is implemented in such a way that simultane- 



ously several solutions corresponding to local and global 
minima of the free energy can be found. The obtained so- 
lutions can then be ranged by their free energy or by the 
concentration profile. Even in the restricted ID geome- 
try of flat phospholipid bilayers with periodic boundary 
condition the method finds several self-assembled struc- 
tures: homogeneous solution of lipids, one single layer, 
formed either at the center or at the edge of the box and 
two or even three layers in the box. Each structure has 
its own free energy which allows to determine the equi- 
librium properties. 

The detailed microscopic information concerning the 
mechanic and thermodynamic properties of the simulated 
flat bilayers can be compared directly against the exper- 
imental data for phospholipid bilayers^—. Thermody- 
namic properties of the membranes obtained from the 
measurements are the thickness of the membrane (deter- 
mined as a distance between the midpoints of the outer 
slopes of the H beads concentration profiles), the full 
thickness of the bilayer region (determined as a thick- 
ness of the region where the total concentration of the 
lipid molecules is larger than zero) , the thickness of the 
hydrophobic core (measured as a distance between the 
midpoints of the slopes of the T beads concentration 
profiles), the distances between different groups of the 
phospholipid molecule, membrane rigidity measured by 
the compressibility modulus and the interfacial area per 
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FIG. 3. Detailed equilibrium concentration profile of the 
phospholipid bilayer obtained with the 44-beads model. 

lipid molecule. Thus, the coarse-grained parameters of 
each model of phospholipid molecule are chosen in such 
a way that all these microscopic properties correspond 
to the experimental values of the DMPC phospholipid 
bilayers. The parameters of these three models are sum- 
marized in Table HI 

The most detailed 44-beads model provides all infor- 
mation about the composition of the membrane as well 
as the position of different groups of the molecule in the 
bilayer. The concentration profile of the bilayer at equi- 
librium is shown in Figure [3l The profiles are normalized 
to unity by division by the occupied volume fraction <j)Q. 
One can see, that the terminal groups of the lipid tails are 
situated in the center of the layer, which is surrounded 
by the area filled with other carbon groups of the tails. 
The groups of hydrophilic units representing the "heads" 
arc in the surface layer. Such position of the groups is in 
agreement with the experimental data for such system^. 
It is noteworthy, that the concentration profile of choline 
group has two peaks reflecting the position of phospho- 
lipid heads in the layer: approximately half of the heads 
are oriented perpendicular to the layer plane, while the 
other half is parallel to the layer plane. This is the out- 
come of the model of the phospholipid molecule and the 
parametrization of the interactions. This can certainly be 
improved once the reliable experimental data or results 
of MD simulations concerning the equilibrium position of 
the heads in the layer is available. 

One of the advantages of the SCMF method is the ac- 
curate measurements of the free energies of equilibrium 
structures. In our model this corresponds to the free en- 
ergy of formation of flat phospholipid layer. It does not 
include the entropy of thermal undulations in the perpen- 
dicular direction. In particular, undulations may result 
in the broadening of the phospholipid membrane. How- 
ever, they become important on the scales much larger 
than the size of the simulation box and, in principle, can 
be included as an additional contribution to the free cn- 



FIG. 4. Free energy per lipid / vs. interfacial area per 
lipid A for the 44-beads lipid model, calculated from solu- 
tions with one (red points) and two (green points) bilay- 
ers in the simulation box. Parameters of the fitting curve 
y — ao + aix + a2X 2 are do = — 61.3±0.5, a\ = — 0.164±0.010, 
a 2 = (1.17 ±0.04) x 10~ 3 . 

ergy of the layer. 

The rigidity of the membrane can be obtained from the 
free energy of the bilayer. Changing the number of lipids 
in the box, we get the free energy per lipid molecule / as 
a function of the interfacial area per lipid A (Figure |4} . 
Although the error bars are quite large because of the 
large conformational space with respect to the sampling 
used for the calculations, the energy points can be quite 
good approximated with the polynomial of the second or- 
der. This curve allows us to calculate the compressibility 
of the layer with respect to compression and stretching 
in the lateral direction. The second derivative is related 
to the compressibility modulus defined in cq. (fl~9|) . The 
resulting compressibility modulus, the thickness of the 
membrane, positions of various groups and the area per 
lipid are shown in Table Ql] The obtained values agree 
quite well with the experimental data. However, if better 
agreement is needed, the parameters of the models can 
be tuned by the additional series of run of calculations 
with a large number of conformations in the sampling. 

The free energy per lipid / (Figure U) has a minimum 
at A « 60 A 2 , which corresponds to the experimental 
value of the interfacial area per lipid in the bilayer in 
the fluid phase. Howether, the values of A below 50A 2 
cannot be obtained with the same set of parameters, be- 
cause this range corresponds to dense packing in the 
core of the layer. Further decrease of interfacial area 
per lipid is accompanied with structural rearrangements 
of the phospholipid environment and strong correlations 
between the neighboring molecules. These correlations 
induce phase transitions in lipid membranes^, the main 
transition to the gel phase occurs ati« 40A 2 . However, 
such correlations are well beyond the limit of validity of 
the mean field theory, where the movements of individual 
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TABLE II. Comparison of the equilibrium properties of phospholipid bilayer obtained with three lipid models with experimental 
data and full atomistic MD simulations of the DMPC lipid bilayer. 





44B 


10B 


3B 


DMPC 


Membrane thickness (A) 


35 


36 


45 


44.2 a 


Full thickness of the bilayer region (A) 


51 


44 


56 


53 6 


Thickness of hydrophobic core (A) 


24 


24 


28 


26. 2 a 


Distance between heads (A) 


30 


28 


36 


36" 


Distance between phosphate groups (A) 


30 






32 6 


Distance between choline groups (1st peak) (A) 


22 






37" 


Distance between choline groups (2nd peak) (A) 


42 






37 M 


Distance between glycol groups (A) 


26 






27 6 


Distance between COO groups (A) 


24 






25 6 


Thickness of terminal groups region (A) 


6 






10 b 


Interfacial area per lipid (A 2 ) 


70 ±5 


65 ± 11 


60 ±2 


59.6 a 


Compress, constant (dyn/cm) 


135 ± 11 


560 ± 120 


269 ± 11 


257 c 



a Experimental data by Nagle and Tristram- NagleP ~ 
b Calculated from MD simulations d ata bvlDamodaran and Mer3 ^ 
c Experimental data bv lMathai et al^ 
d Average distance between choline groups. 




FIG. 5. Free energy per lipid / vs. interfacial area per lipid A 
for the 10-beads lipid model, calculated from solutions with 
one (red points), two (green points) and three (blue points) 
bilayers in the simulation box. The inset represents the en- 
larged region of the minimum. Parameters of the fitting curve 
y — a + aix + a 2 x 2 are ao — —3.0 ± 3.0, <zi = —0.68 ± 0.08, 
a 2 = (5.2 ±0.6) x 10~ 3 . 

molecules are supposed to be uncorrelatcd ^j. 

Similar calculations carried out for the 10-beads model 
show better reproducibility due to considerably reduced 
conformational space. In fact, 1 million of conformations 
is enough for calculations with different sampling to al- 
most coincide in one curve with very small error bars 
(Figure [5J . The absolute values of the free energy / are 
shifted with respect to that of 44-units model, which is 
the consequence of different sizes of the units and the in- 
teractions parameters. However, the energy differences, 



FIG. 6. Free energy per lipid / vs. interfacial area per lipid A 
for the 3-beads lipid model, calculated from solutions with one 
(red points) and two (green points) bilayers in the simulation 
box. The inset represents the enlarged region of the minimum. 
Parameters of the fitting curve y — ao + a\x + a 2 x 2 are ao = 
-1.4 ± 0.2, ai = -0.325 ± 0.008, a 2 = (2.72 ± 0.06) x 10~ 3 . 

the density profiles and the compressibility modulus are 
very close to that obtained with the previous model (Ta- 
ble [IlJ) . It is interesting to compare the energies of one 
and two layers in the box, red and green curves in Fig- 
ure [SJ correspondingly. The energies are quite similar at 
low densities of the layer and diverge significantly at high 
densities. This is consequence of the compression of the 
two layers due to lack of space in the box for two layers of 
equilibrium size. The energy per lipid corresponding to 
the maximal compression of two layers in the box coin- 
cides with the energy of a compressed layer. This fact re- 
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FIG. 7. Mean field equilibrium "snapshots" (set of most-probable conformations of the lipids) and the concentration profiles of 
tails (T, red lines), heads (H, blue lines) and the phospholipids (T+H, black lines) obtained with 44-, 10- and 3-beads models. 



fleets the consistency of our picture. On the other hand, 
high compression of a single layer, A below 50 A 2 , results 
in gradual splitting of the layer in two with smooth re- 
arrangement of the lipid molecules. Some phospholipids 
flip and orient their "heads" in the center of the mem- 
brane. Thus, a single bilayer smoothly become broader 
and change its appearance to remind something between 
a single bilayer with a layer of heads in the center and 
two tightly compressed bilayers. This process is reflected 
in the flat part of the energy per lipid curve for A below 
50 A 2 , where the distinction between one layer and two 
layers in the box can hardly be done. 

The energy curve for the 3-beads model is presented in 
Figure [6] This is the simplest phospholipid model with 
only 3 beads, thus the conformational space is very lim- 
ited and the results obtained with the sampling of 1 mil- 
lion of conformations represented here are close to these 
obtained with small sampling about 100-300 thousands. 
The difference between the energies of one and two lay- 
ers in the box, noted for 10-beads model, is even more 
pronounced. 

The mean field " snapshots" of the most probable con- 
formations of the layer in equilibrium and the corre- 
sponding density profiles for three models are presented 
in Figure [Jj It shows how the coarse-graining and the 
decreasing of the model details influence the layer ap- 
pearance. However, despite obvious differences in details 
caused by the differences in the molecule structures, the 
profiles looks quite similar, especially for the most coarse- 



grained 10- and 3-beads models. In case of the 44-units 
model, there is a small fraction of the hydrophobic T 
units at the surface, reflecting the structure of the phos- 
pholipid with some hydrophobic T groups in the head. 
Furthermore, different aspect ratio between T and H 
units in the models (4 : 1 in 10-beads and 2 : 1 in the 3- 
bcads model) lead to the corresponding difference in the 
relative density of T and H units. This correspondence 
is not valid for the 44-beads model, where the beads do 
overlap and some parts of the molecule are hidden from 
the solvent. The ratio between the areas under the T 
and H curves is 4.7 : 1 which is different from the ra- 
tio between numbers of T and H beads in the molecule 
(3.4:1). 

Finally few words should be said about the speed per- 
formance of three models. We run six series of calcula- 
tions with different conformational sampling of 1-3 mil- 
lions of conformations. The box is divided in ~ 100 par- 
allel layers and the number of lipids in the box is varying 
from 300 to 700 with the step 10. The program runs in 
parallel on 32-cores AMD machine with 128Gb of RAM. 
The most time consuming process is the solution of equa- 
tions while the generation of the sampling at the begin- 
ning of calculation is relatively fast and almost does not 
affect the speed of calculation. In turn, the solution of 
the equations depends on the number of conformations 
in the sampling, the number of cells and the number of 
attempts to solve equations with different initial condi- 
tions. Since three models differ in the number of beads, 
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implying different conformational space, the same accu- 
racy of calculations in different models is achieved with 
different number of conformations in the sampling. The 
44-beads model require several millions of conformations 
and 48 hours of calculations, the 10-beads model require 
one million of conformations and 17 hours, while the 3- 
beads model can get reliable results with a hundred thou- 
sands of conformations in few hours. Thus, the simplest 
3-beads model is the fastest model which can reproduce 
the essential equilibrium properties of the phospholipid 
bilayer. 

V. CONCLUSIONS 

In this work we report a computational tool provid- 
ing relatively fast and stable solution of the equations 
of the SCMF theory in different geometries and different 
molecule structures. The SCMF theory is applied for the 
first time to model the DMPC phospholipid membranes 
with realistic thermodynamic properties. We show that 
all three coarse-grained models of phospholipids with the 
present parameters can adequately describe the equilib- 
rium properties and the molecular structure of DMPC 
phospholipid bilayers in fluid phase. Among the essential 
properties of the phospholipid bilayers that correspond to 
the experimental data are the thickness of the membrane, 
the positions of different groups of the lipid in the bilayer, 
the compressibility of the membrane and the equilibrium 
interfacial area per lipid. The most detailed 44-beads 
model can reproduce the structure of the phospholipid 
layer with maximum details. The resulting density pro- 
files of the groups of the phospholipid molecule can be 
used for the structural modeling in experimental tech- 
niques that require the molecular structure of the layer. 

However, if one is interested in the average proper- 
ties of the lipid bilayer, the simplest 3-beads model can 
be used to model the thermodynamics and the essential 
properties of self-assembled phospholipid bilayers. This 
model is faster and computationally more attractive than 
44- and 10-beads models. Furthermore, we suggest that 
the 3-beads model with the parameters found in this work 
is used for simulation of more complex bilayer structures. 

ACKNOWLEDGMENTS 

The authors are grateful to Josep Bonet and Allan 
Mackic for introducing into the method SCMF the- 
ory, Nigel Slater for his hospitality during the visit of 
the University of Cambridge and stimulating discus- 
sions. The authors gratefully acknowledges financial help 
from Spanish Ministry of education MICINN via project 
CTQ2008-06469/PPQ and the UK Royal Society for the 
International Joint Project. 



Appendix A: Computational Details 

Here we discuss the implementation details of the al- 
gorithm underlying the computational realization of the 
SCMF theory. The off-lattice realization of the program 
allows for conformations of the molecules in a real space. 
The algorithm consist of three independent parts: the 
generation of the sampling of conformations of molecules 
of any structure, the discretization of space according 
to the geometry of the problem and the solution of the 
SCMF equations (fT2")l . The independence of the parts in- 
sures the universality of the method which can be applied 
for molecules of any structure and their self-assembly in 
objects of any geometry. 

1. Conformational sampling generation 

First step is the generation of a representative set of 
conformations of a single molecule. This can be done 
with different techniques, including bead-to-bead chain 
growth or MC simulations. The conformations of a single 
molecule are generated once before solving the equations. 
They are stored in the RAM memory and not changed 
during the solution of equations, while the probabilities of 
each conformation are recalculated with the fields. This 
static memory allows for highly efficient parallelization of 
the program within the shared memory OpcnMP plat- 
form, since the processors communicate very little be- 
tween each other. The generation of conformations in 
the sampling is realized as follows: 

1. Generation of a new conformation r Q of the 
molecule of type a with random position of the 
first bead inside the box. We set initial values of 
the conformation statistical weight w a (T a ) = 1 and 
the intramolecular energy H™ tra (T a ) = 0. 

2. Iterative addition of new beads to the growing 
molecule and joining them one by one subject to 
self-avoidance condition. Depending on the final 
objective, one can use: 

(a) Self-avoiding growth with the MC equilibra- 
tion. The chain is generated as a self-avoiding 
off-lattice walk and the energy of interaction 
of the newly joined beads with the rest of 
the chain is accumulated in H™ tra (T a ). The 
resulting chain is then equilibrated with MC 
crankshaft moves and reptations until the full 
equilibration of the configuration. In this case 
all conformations have the same statistical 
weight and w a (T a ) = 1. 

(b) Rosenbluth chain growth^. A new bead is 
placed at a fixed distance from the previous 
bead either at random angle (freely joined 
chain) or with additional restrictions on the 
angles between the beads (chain with angle 
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restrictions). Rosenbluth algorithm generates 
conformations with biased probability distri- 
bution. This bias is then calculated in the 
conformation weight and removed when the 
probabilities of each conformation are cal- 
culated. In order to calculate the Rosen- 
bluth weight, one should try N tr iai times to 
place the new bead at a fixed distance from 
the previous bead and calculate the number 
Naiiowed of successful positions, allowed by 
self-avoidance condition. If N a iiowed > 0, 
a new position is accepted with the weight 
]-/N a uowed, the energy of interaction is accu- 
mulated in H l ™ tra (T a ) and the weight of the 
conformation w a (T a ) is multiplied by the fac- 
tor 1/N allowed- If there is no possibility to 
place the bead, N a u owec i = 0, the generation 
restarts from the beginning. The Rosenbluth 
algorithm is simple and efficient enough for 
generation of phospholipid molecules. 



2. Discretization of space 

Since the solution of equations of the SCMF theory in 
the integral form (fT2|) . (|13[) is not practical, the equations 
are simplified by the replacement of the integrals with the 
sums. The space is discretized according to the geome- 
try of the system and the symmetry considerations. The 
integrals over the spatial coordinates, dr, are replaced by 
the sums over auxiliary cells i with the same values of 

the mean fields: /c^(r)\ = and A(r) = Aj. The dis- 
cretization leads to the following approximate relations 

J <(r Q ,r) (4(r)) dr*J2 c h J<(T a ,r)dr 
J X(r)(j) a (r a ,r)dr w J (j) a {T a ,r)dr 



(Al) 



where ^ i denotes the sum over all auxiliary cells, while J. 
is the integral over the i-th cell. Because of the structure 
of the equations, the integral over the single i-th cell can 
be evaluated once in the beginning of calculations, while 
the sums (|A1[) are calculated during the solution of the 
equations. Although the number and the geometry of 
the cells can be arbitrary, their choice will restrict the 
geometry and the resolution of the resulting solutions. 
Thus, it is important to choose them accurately for a 
particular problem. In case of flat bilayers we choose a 
planar geometry with the cells parallel to the layer, thus 
we can only obtain the solutions that can lie in the layer. 

The integrals over the conformational space of single 
molecules are replaced by the sums over the finite sam- 



plings generated in the first step 

4>ai = ^ </>m(r a )p a (r a ), 



(A2) 



where the mean field volume fraction <f> a i and the mean 
field concentration c a ai of a-beads of the molecule of type 
a in the cell i are related to the corresponding volume 
fraction (j)ai(Xa) and the concentration c^iTa) of the 
conformation r Q . 

The integral expression (fTTj) for the p a (T a ) yields in 
form 



Pa(r„) 



1 



Z a W a (T a ) 



exp(- J ff- t ™(r Q )- 



M 



(A3) 



(r Q ) c ^ + 2^A^ m (r Q ) 



where w a (T a ) is the Rosenbluth weight calculated during 
the generation of the chains, and the constant Z a insures 
the normalization of the probability p Q (r Q ) = 1. 

The expression for probabilities (|A3|) is accompanied 
with the following relations 

M 

v s X l « ln^ + ^2(Nfi~ 5 sf3 ) 

0=1 a 
M 

<t>si = 00 - ^ N a (j) m 



a— l,a^s 



(A4) 



&(Ta)= fu a a (T a ,r)dr, 

J i 

v a i(T a ) = Vi4> a i(r a ) = I Q (r a ,r)dr 



The quantities v ai (T a ), e^(r a ) and c^(r a ) are calcu- 
lated using the MC integration technique. The volumes 
v a i(Ta) are calculated as follows: one throws randomly 
Ntriais points inside the volume occupied by the bead 
and determine the cell where each of these points hits. 
Each point contributes to the volume in the cell with the 
value v a j (Ntriais Ni„t), where Ni n t is the number of the 
beads hit by the point and v a is the volume of the bead. 
The volume fractions (j> a i(T a ) are related to the volumes 
Vai(T a ) as (j) a i(T a ) = v a i(T a )/Vi, where Vi is the volume 
of the i-th auxiliary cell. The concentrations c^fTo.) 
are calculated similarly, except that each random point 
contributes to c^(T Q ) with the value 1/ Ntriais- The in- 
teraction energies e°j(r a ) arc calculated with the help 
of the Ntriais points thrown for every a at a distance r 
from the center of the bead b such that r m ;„ < r < r max , 
where r m in is the minimum distance at which the bead of 
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type a can approach to the b bead, and r max is the max- 
imum distance at which the interaction potential acting 
between a and b beads is non-zero. If the bead of type 
a in the position of the trial point does not have any in- 
tersection with the beads of the chain, its contribution 
to the e a ai (T a ) is §tt (r 3 max - r s min ) U ab (r)/N trialf , where 
i is the number of the auxiliary cell where the trial point 
hits, and U ab (r) is the pair potential acting between the 
beads of types a and b. Since the solvent is implicit in 
our model, the quantities describing the interaction 
between the solvent and the bead of type a have to be 
treated separately. We evaluate them by a simple ap- 
proximate expression 

<%i « ^O-L - (r s + r a ) 3 )e sa (A5) 

where we assume the square well potentials between the 
solvent molecules s and the beads a with energy e sa at 
the distance closer than r sa . Here the radius of solvent 
is r s and the radius of the a bead is r a . 

It is noteworthy that the incorporation of the Roscn- 
bluth weight w a (T a ) into the expression for probabilities 
(|A3|) may lead to confusion. Although most of the av- 
erages can be calculated with the expressions similar to 
(|A2p . the averages of expressions containing explicitly the 
probability p a should be treated separately. For exam- 
ple, the entropy term in the free energy (|10j) leads to the 
additional term when we replace the integral with the 
sum 

N a (in *f*A a ) = N a E rct Pa (T a ) In p ° (r ; )jV ° A a + 

N »J2r a Pa(T a )\iiw a {r a )x a (A6) 

where x a is the number of conformations in the sampling 
of the molecule of type a. 

The expressions (|A2|) . (fA3j) and (|A4|) form a closed 
set of algebraical equations with respect to variables 4> a i 
and c a ai which can be solved numerically. We must admit 
however, that the solution of these non-linear equations 
is not a straightforward task. 

3. Numerical implemetation 

In this work we have implemented the SCMF theory 
in C++. The program is based on Newton-Rapson tech- 
nique which provides the efficient and stable solution of 
the SCMF equations. It has a modular structure which 
makes it very flexible for modification, further develop- 
ment and extension for new applications. 

The first module is the sampling generator, which gen- 
erates the molecules conformations and calculates the 
volumes and interactions according to the Appendix IA II 
The code for generation of molecules with different struc- 
tures can be added to this module or changed indepen- 
dently from other parts of the program. 

The second module specify the geometry of the prob- 
lem: ID spherical, ID planar, 2D planar, 2D cylindrical, 



or 3D cubic. It contains the description of the geometry 
and divide the simulation box into set of auxiliary cells. 
The modification of this part of code provides an easy 
way to implement different types of the system geome- 
try. It is worth to note here one technical improvement 
that significantly increase the performance of the calcu- 
lations. Since the molecules have small size compared 
to the size of the box, their contribution to most of the 
auxiliary cells is zero, thus matrices e^(r Q ), c^j(r Q ) and 
v a i(Ta) have many zero values. In this second module of 
the program the simple compression of these data is re- 
alized via elimination of zeros from these matrices. This 
helps to save computer memory and, on the next step, 
speed up the calculations by elimination of the sums over 
the cells containing zeros. 

The core of the third module is the modified Newton- 
Rapson solver for the solution of SCMF equations. It in- 
volves the iterative solver which calculates on each step 
the matrix of analytical derivatives of equations IA2I with 
respect to unknown variables along with the analytical 
expressions including the sums over the conformational 
samplings. The next approximation of the solution is de- 
termined by the inversion of this matrix until the required 
accuracy is reached. To avoid occasional instability of 
the Newton-Rapson scheme, the algorithm artificially re- 
strains too big steps which may drop the next approxi- 
mation outside a reasonable range. The algorithm is able 
to find several solutions simultaneously by repeating the 
calculations several times with different initial conditions. 
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